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We propose a general strategy to develop accurate Force Fields (FF) for metallic systems derived 
from ah initio quantum mechanical (QM) calculations; we illustrate this approach for tantalum. 
As input data to the FF we use the linearized augmented plane wave method (LAPW) with the 
generalized gradient approximation (GGA) to calculate: 

(i) the zero temperature equation of state (EOS) of Ta for bcc, fee, and hep crystal structures for 
pressures up to 500 GPa. 

(ii) Elastic constants. 

(iii) We use a mixed-basis pseudopotential code to calculate volume relaxed vacancy formation 
energy also as a function of pressure. 

In developing the Ta FF we also use previous QM calculations of: 

(iv) the equation of state for the A15 structure. 

(v) the surface energy bcc (100). 

(vi) energetics for shear twinning of the bcc crystal. 

We find that with appropriate parameters an embedded atom model force field (denoted as qEAM 
FF) is able to reproduce all this QM data. Thus, the same FF describes with good accuracy the bcc, 
fee, hep and A15 phases of Ta for pressures from ~ — 10 GPa to ~ 500 GPa, while also describing the 
vacancy, surface energy, and shear transformations. The ability of this single FF to describe such a 
range of systems with a variety of coordinations suggests that it would be accurate for describing 
defects such as dislocations, grain boundaries, etc. 

We illustrate the use of the qEAM FF with molecular dynamics to calculate such finite temperature 
properties as the melting curve up to 300 GPa; we obtain a zero pressure melting temperature of 
Tmeit = 3150 ± 50 K in good agreement with experiment (3213 — 3287 K). We also report on the 
thermal expansion of Ta in a wide temperature range; our calculated thermal expansivity agrees 
well with experimental data. 
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I. INTRODUCTION 

Despite decades of experimental and theoretical re- 
search on the mechanical properties of materials many 
questions remain open, particularly in the relation be- 
tween atomistic processes (involving dislocations, grain 
boundaries, cracks, their mobility and interactions 
thereof) and the macroscopic behavior (plastic deforma- 
tion, failure, etc.). Macroscopic plasticity and failure 
are well characterized experimentally and described using 
a variety of mesoscale and macroscale models with pa- 
rameters obtained empirically. These models and their 
parameters should ultimately be derivable in terms of 
the fundamental physics of atomic interactions as de- 
scribed by quantum mechanics (QM). Unfortunately de- 
spite the enormous progress in ab initio QM, such cal- 
culations are too computationally demanding to study 
directly most processes relevant to plasticity and failure. 



To do so would require that millions of atoms be de- 
scribed for nanoseconds or longer, calculations that are 
impossible to consider today. In order to bridge this gap 
between atomic interactions and the mechanical proper- 
ties of macroscopic systems we use first principles QM 
data to derive a FF with which energies and forces can 
be calculated in a computationally efficient way given 
only the atomic positions. This allows us to use clas- 
sical molecular dynamics (MD) to simulate the various 
atomistic processes governing the mechanical and thcr- 
modynamical properties of materials. Since the FF is 
derived by fitting a wide range of QM data, we expect to 
obtain an accurate description of the atomic interactions 
which, with molecular dynamics, can provide insight and 
constitutive equations to be used in the mesoscale and 
macroscale models of plasticity and failure. We illus- 
trate this procedure here for Ta. Ta is chosen because it 
exhibits only a single crystalline phase over the interest- 
ing range of temperatures and pressures. This makes the 
validation of theoretical predictions against experimental 
data under a wide range of conditions more clear. 
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Several authors have used ab initio methods to study 
various mechanical properties of Ta. Soderlind and Mo- 
riarty used the full potential linear muffin-tin orbital 
method within the GGA approximation and with spin 
orbit interactions to calculate different zero temperature 
properties of of Ta, including the EOS of different [Crys- 
talline phases, elastic constants, and shear strength.^! Va- 
cancy formation migration|^ergies have also been calcu- 
lated-from first principles Elj Recently, Ismail-Beigi and 
Ariasu calculated ab initio dislocation core energy and 
structures for Ta and Mo. 

On the other hand, several FF have been developed for 
bcc metals. The first principles based multi-ion analytic 
model generalized pseudopotential (MGPT) force field 
has been used to calculate several mechanieaL and thcr- 
modynamical properties of various metals,El Q including 
Ta. One of the most popular many-body force fields for 
metals is the Embedded Atcup Model (EAM), proposed 
by Daw and Baskes in 1984.B Most EAM FF have been 
based on experimental data regarding structures at and 
near equilibrium [for bcc metals, see for example Ref. 
|-|n|. One of the main advantages of EAM force fields 
is that they are computationally very efficient which al- 
lows MD simulations of large systems for long times. 
The modified embedded atom model (MEAM), which in- 
cludes angular dependence of the eloetronic density, was 
developed by Baskes and coworkerst3 and appliedJo a 
variety of materials including bcc transition metalaij. 

In this paper we propose a strategy to derive accurate 
many body FF based on QM calculations; these FF can 
be used to simulate the dynamics of systems with millions 
of atoms. We use accurate QM (LAPW GGA method) 
to calculate various mechanical properties of Ta which 
require a small number of atoms: namely the zero tem- 
perature EOS for bcc, fee and hep crystalline phases, and 
elastic constants. We also use a mixed-basis pseudopo- 
tential method to calculate vacancy formation energies. 
We find that we can describe all this first principles data 
using a classical many body EAM FF (named qEAM 
FF) with good overall accuracy. We then illustrate the 
use of the qEAM FF with MD simulations to study vari- 
ous properties as a function of pressure and temperature; 
such as the melting curve of Ta up to pressures of ~ 300 
GPa, and thermal expansivity. 

This paper is organized as follows. In Section II we 
present first principles results for the equation of state 
for bcc, fee and hep phases, elastic constants and volume 
relaxed vacancy formation energy and enthalpy for a wide 
pressure range. In Section III we develop the qEAM FF 
based on ab initio results. In Section IV we calculate the 
thermal expansivity of Ta using the qEAM FF and ab 
initio calculations. Section V presents the calculation of 
the melting curve of Ta for pressures up to 300 GPa using 
molecular dynamics. Finally, in section VI, conclusions 
are drawn. 



II. QUANTUM MECHANICS RESULTS 

We computed the static equation of state of Ta for 
different crystalline phases using the-, [linearized aug- 
mented plane wave (LAPW) methodo^E^I LAPW is an 
all-electron method, with no essential shape approxima- 
tions for the charge density or potential, and is easily 
converged. The bp, 4/, 5d and 6s states were treated as 
band states, and the deeper states were treated as soft 
core electrons. Here we used the Perdew, Burke, and 
Ernzerhof (PBE) implementation of the generalized gra- 
dient approximationE3 for the exchange-correlation po- 
tential. A 16x16x16 special k-point meshliJ was used, 
giving 140 k-points within the irreducible Brillouin zone 
of the bcc lattice. Tests demonstrated convergence with 
this mesh. The convergence parameter KK^max was 9 giv- 
ing about 1800 planewaves and 200 basis functions per 
atom at zero pressure. 

Total energies were computed for bcc, fee and hep 
phases at 20 volumes. For the fee and hep phases, 
12x12x12 and 16x16x12 k-point meshes were used for 
Brillouin zone integrations giving 182 and 180 k-points 
within the irreducible zone respectively. For the hep 
phase, the ideal c/a ratio was used, and at two differ- 
ent volumes also c/a was optimized. We found that 
the change in the energy due to the this optimization 
is less than 40 meV/atom around zero pressure and de- 
creases by pressure. We have tabulated the ab initio 
energy-volume data for the different phases and have 
made them available as supplementary material, see Ref. 
|l^ . We have fitted our energy-volume data to Rose's 
universal equation of statejl3 the obtained zero pressure 
volume (Vq), zero temperature bulk modulus (Bt), and 
its derivative with respect to pressure (Brp) are shown in 
Table ||. We also show in Table | the results obtained by 
Soderlind and Moriarty using full potential linear muffin- 
tin orbital method within the GGA approximation and 
with spin orbit interactions (denoted as FP LMTO GGA 
SC), and[«)om temperature experimental values by Cynn 
and YooE3 Our LAPW calculations of the bcc equation 
of state agree well with the experimental values and pre- 
vious theoretical calculations. 

Static elastic constants [Cs=(cii-Ci2)/2 and C44] were 
obtained from strain energies by straining the bcc cell 
with volume conserving tetragonal and orthorhombic. 
We calculate Cs using tetragonal strain of the cubic bcc 
lattice: 

a =a(l 4-6,0,0), 
b = a(0,H-e,0), 

c = a(0,G,l/(l + e)2), (1) 

where a is the cubic lattice constant of the system, a, b, 
and c are the lattice vectors and e is the strain. Cs is 
related to the quadratic term of the strain energy: 

E{e) = E{e = 0) + &V{e = 0)c,e2 + Qie'), (2) 
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where E{e = 0) is the energy of the unstrained system 
and V{e = 0) is its volume. Similarly C44 is obtained 
from the orthorhombic strain: 



a = a(l,e,0), 
b = a(e,l,0), 
c = a(0,0,l/(l-e2)). 



(3) 



the shear constant C44 is obtained from: 

E{e) = E{e = 0) + 2V{e = 0)0446^ + 0{e^). (4) 

The convergence of strain energies with respect to the 
Brillouin zone integration was carefully checked; we used 
16x16x16 k-points meshes in the full Brillouin zone giv- 
ing 344 and 612 k-points within the irreducible Brillouin 
zone of tetragonal and orthorhombic lattice respectively. 
Cs and Bt = {cu + 2ci2)/3 were used to calculate cn 
and C12; the resulting zero pressure and zero tempera- 
ture elastic constants are shown in Table |. The elastic 
constants Cs and C44 as functicps of pressure are avail- 
able as supplementary materialllS, see also Ref. The 
zero pressure values and initial slopes are in good agree- 
ment with the experimental data of Katahara et al.E3 We 
find that C44 shows a change in behavior at ~150 GPaEiJ 
which is probably due to the eleetronic transition also 
evident in the equation of state.l2a The band structure 
and den^ty of state show a major reconfiguration with 
pressure.E3 Our results indicate that the elastic constants 
can be much more sensitive to changes in the occupied 
states below the Fermi surface than the equation of state, 
where changes were much more subtle. j— , 

The vacancy formation energy was determinedHJ us- 
ing supercells of 16 and 54 atoms and a mixed- 
basis pseudopotmitial (MB-PS) codeBj Using the GGA 
approximation ,0 the zero pressure volume relaxed va- 
cancy formation energy was determined to be 3.25 eV 
for a 54 atom supercell, and 3.26 eV for a 16 atom su- 
percell, indicating convergence. The effects of structural 
relaxation on vacancy formation energy are discussed in 
Ref. A comparison with recent ab initio and experi- 
mental results will be presented in the next section. The 
ab initio data used to calculate vacancy formation ener- 
gies is also available as supplementary material. t3 



III. qEAM FORCE FIELD 

In order to predict mechanical properties of materials 
and processes like phase diagrams, dislocations structures 
and mobility, mechanical failure, etc. it is important to 
have accurate classical force fields to describe the atomic 
interactions. With MD simulations it is then possible to 
study large systems (millions of atoms) for relatively long 
times (ns). 

One of the most popular many-body force fields for 
metals is the EAM, proposed by Daw and Baskes in 
1984.eI This approach is computationally efficient and has 



been used successfully for numerous applications, like cal- 
culation of thermodynamic functions, liquid metals, de- 
fects, grain boundary structure, fracture, etc, see for ex- 
ample Ref. |ll|,||. 

The EAM implementation that we have used is based 
on the one proposed by Chantasiriwan and Milstein.EI 
This from of EAM was chosen becauiie it can describe 
third order elastic constants correctlyQ thus being useful 
in a wide range of strains. The energy of a given atomic 
configuration with atom positions {vi} is given by: 



N 



U{n} = J2F{p.,) + J2<l>{nj), 



with 



(5) 



(6) 



where F{p) is the embedding energy, pi is the total 
"electronic density" on the atomic site i, /(r) is the 
electron density function, 4>{r) is a two-body term and 

Tij — I 

Following Ref. || we took the two-body term as: 



<t>{r) 







otherwise. 



(7) 



the factor (r — r™)^ ensures that the two body term and 
its first three derivatives with respect to r vanish at the 
the cut-off radius (r^). The optimized form of the two 
body term shows short range repulsion and longer range 
attraction. 

The "electron density" is: 



f{r) 



1 + ai cos(-rn73 ) 02 sin(-r^) 



(8) 



where V is the volume per atom of the system. The 
parameters of the "electron density" are volume depen- 
dent, but structure independent. The importance of the 
oscillatory behavior of the "electron density function" in 
embedded atom model-like force fields is related to their 
ability to correctly account anharmonicities.u 

Finally the embedding energy as a function of the 
electroniCp, density is obtained from the reference bcc 
structure:la 



F{p)^URose{V)~Y.' 



(9) 



where the sum is made for the perfect bcc striicture and 
URoseiy) is Rose's universal equation of state:Ej 

URoseiV) = ~E,oh{l + a* + /3a*3 + Ua*^)e-''' (10) 

where a* = {a - ao)/aoA and A = (£^co?i/9V"oi?T)^^^- 

We we optimized the parameters in the qEAM FF us- 
ing the following data: 
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(a) Zero temperature energy-volume and pressure- 
volume curves for different crystal structures (bee, fee 
and A15) in a wide pressure range, from — 10 GPa to 
~ 500 GPa. For the bee and fee structures the data from 
Section II was used, while for A15 pye used the results 
obtained by Soderlind and MoriartyJd 

(b) Zero temperature, zero pressure elastic coefficients, 
shown in section II. 

(c) Vacancy formation energy at zero pressure. 

(d) Energetic of homogeneously sheared bcc crystal, 
from Ref. [|. 

(c) Unrelaxed (100) surface energy from Ref. 

We fit the parameters entering the qEAM FF energy 
expression to the training set using an optimization algo- 
rithm based on simulated annealing. We define an error 
function of the form: 



(Q 



qEAM 



-\target\2 



{Q 



target\2 



(11) 



where the sum runs over the different quantities to be fit- 
ted, Ql^^'^ is the value given by the qEAM FF for quan- 
tity i, Q*°-^9et ^j^g target quantity and Ci are weight 
factors. The EOS information of the different phases ac- 
count for ~ 78% of the total cost; elastic constants and 
vacancy formation energy account for ~ 9% each one; 
surface energy represents ^ 3% of the total cost and en- 
ergy of sheared bcc crystal represent ~ 1%. 

In Tables and ^ we show the optimized qEAM pa- 
rameters. 

A detailed comparison between the qEAM FF and the 
data it was fitted to is shown in the following subsec- 
tions, together with other ab initio and experimental re- 
sults. We also show the calculation of related quantities 
obtained using the qEAM FF via MD simulations. 



A. Equation of state and elastic constants 

The most important quantities used to develop the 
qEAM FF are zero temperature EOS for different crystal 
phases of Ta in a wide pressure range. We used energy- 
volume and pressure-volume for bcc Ta from Section II; 
fcc-bcc energy difference and fee pressure for different 
volumes from Section II; and first principles Al5-bcc en- 
ergy difference obtained by Soderlind and Moriartyfil us- 
ing full potential linear muffin-tin orbital method within 
the GGA approximation with spin orbit interactions. Ta 
is a bcc metal and no pressure-induced phase transition 
to other solid structure has been found experimentally or 
theoretically. Nevertheless, using QM it is possible to cal- 
culate the equation of state for different crystalline struc- 
tures although they may not be thermodynamically sta- 
ble. Including the EOS of different phases in the FF de- 
velopment leads to an accurate description of the atomic 
interactions even when the environment of an atom is 
not that of the stable phase; this could play a key role to 
correctly describe defects and non-equilibrium processes. 



In Figure 1 we show energy [Figure 1(a)] and pres- 
sure [Figure 1(b)] as a function of volume for bcc Ta at 
T = K. The circles denote LAPW GGA resuhs and 
the lines the qEAM FF. Figures 2 and 3 show the same 
results for fee and A15 Ta respectively. In Figure 3(a) 
open circles slapw the A15 energy calculated by Soderlind 
and Moriarty,El and the filled circles denote the sum of the 
A15-bcc energy difference from Ref. |l| and our bcc energy 
from section II, these are the quantities the FF was actu- 
ally fitted to. The insets on Figures 2(a) and 3(a) show 
the fcc-bcc and A15-bcc energy difference as a function 
of volume obtained with the qEAM FF (lines) and from 
QM (circles). In Figure 4(a) and 4(b) we show energy- 
volume and pressure- volume curves for hep phase respec- 
tively; circles denote denote the LAPW GGA results of 
section II and the lines show the qEAM FF results. Note 
that although the the hep data was not included in the 
training set for the qEAM FF the agreement is very good. 
Figures 1 to 4 show that the qEAM FF reproduces the 
zero temperature EOS for the four different phases very 
well. 

We also included in the FF training set the ab initio 
elastic constants from Section II at zero pressure. Table 

shows bcc Ta EOS parameters [zero pressure volume 
Vb), bulk modulus (Sr), its derivative with respect to 
pressure (B'rp) and the elastic constants] obtained using 
the qEAM FF together with the QM values from Sec- 
tion II and the ones reported in Ref. |l|. Vq, Bt, and B'j, 
were obtained fitting Rose's universal equation of state 
to the energy- volume data shown in Figure 1. The elas- 
tic constants Cg and C44 were calculated with the qEAM 
FF using the tetragonal and orthorhombic strains shown 
above [equations (1) to (4)]. 

In Figure 5 we show the elastic constants [bulk mod- 
ulus Bt, Cs and C44] as a function of pressure obtained 
with the qEAM FF (filled circles and full lines) and the 
LAPW results from Section II. While the agreement in 
Bt is excellent and that for Cs is good, qEAM FF greatly 
underestimates C44 for high pressures; this problem is 
amplified by the possible electronic phase transition that 
leads to a change of behavior of C44 at ~ 150 PGa (see 
Section II and Ref. |l|). 

In order to estimate the accuracy of the the numerical 
calculation of elastic constants we have calculated cn 
independently using a uniaxial strain: 



a = a(l + e,0,0), 
b = a(0, 1,0), 

c = a(0,0,l), (12) 
the elastic constant cn is obtained from: 

E{e) = E{e = 0) + PV{e = 0)e + ^V{e = 0)cne^ + O(e^), 

(13) 

where E{e ~ 0) is the zero strain energy, V{e — 0) 
is the volume at zero strain, and P is the zero strain 
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hydrostatic pressure. For zero pressure we calculate 
cii = 273.67GPa only 0.4 % higher than the one com- 
puted from Cs [equations (1) and (2)] and Bt (calculated 
form the energy- volume data shown in Figure 1). We 
find a similar agreement under compression: for pres- 
sure P=109.54 GPa (corresponding to a volume per atom 
V{e = 0) = 13.36A3) we obtain cn = 819.90 GPa only 
0.8 % larger than the calculated from Cg and Bt- 

We have calculated the T= 300 K EDS using 
isothermal-isobaric (NPT) MD with a HpoveiEj thermo- 
stat and Rahman-Parrinello barostat.E3 In Table | we 
show the T=300 K zero pressure volume, bulk modu- 
lus and its first derivative with respect, to volume; we 
also show recent compressibility dataE3 obtained in a 
diamond-anvil cell at room temperature and-iultrasonic 
measurements of adiabatic elastic constants£3 

Zero temperature EOS data obtained with tte qEAM 
FF is available in the supplementary material.113 



B. Vacancy formation energy 

We used experimental values for vacancy formation en- 
ergy and cohesive energy in the training set used to fit 
the qEAM FF. The experimental value for_±he relaxed 
vacancy formation energy is Eyac = 2.8 eV.EJ From this 
value we estimated the value of the unrelaxed vacancy 
formation energy to be 3 eV. This unrelaxed value was 
used to fit the qEAM FF. The value for the cohesive 
energy used is Ecoh = 8.10 eV.E3 

The volume-relaxed vacancy formation energy is de- 
fined in computer simulations as: 

iV - 1 

e.ac(-P) = e,ac{N ~1,P)~ -^e^taiiN, P), (14) 

where Cxtai {N, P) is the energy of the perfect crystal with 
N atoms at pressure P and e^aciN — 1, P) is the energy 
corresponding to the N-1 atoms system with a vacancy at 
pressure P where the atomic positions are not allowed to 
relax. The relaxed vacancy formation energy is defined 
in the same way but with eyac{N —1,P) being the atom- 
relaxed energy of the system with a vacancy. 

The volume-relaxed vacancy formation energy ob- 
tained using the qEAM FF is 3.10 eV, in very good agree- 
ment with the target value (3 eV) and only slightly lower 
than the ab initio value 3.25eV. The relaxed vacancy for- 
mation energy obtained with the qEAM FF is 2.95 eV. 
Previous work by Korhonen, Puska and Nieminen, us- 
ing a DFT full potential linear muffin-tin orbital method 
with LDA giveg-|3.49 eV for the unrelaxed vacaacy for- 
mation energy;E2l Satta, Willaime and GironcoliO (using 
DFT plane Waves LDA) calculated 3.51 eV and 2.99 eV 
for unrelaxed and relaxed vacancy formation energies; 
recently Soderlind, Yang, Moriarty and Wills (also using 
the full potential linear mufhn-tin orbital method with 
GGA) obtained 3.74 eV for the volume relaxed vacancy 
formation energy and 3.2 for the relaxed vacancy forma- 
tion energy. Table IV summarizes our vacancy formation 



energies together with previous theoretical results and 
experimental data. 

In Figure 6(a) we show the volume-relaxed vacancy 
formation energy (e.vac) as a function of pressure (P); 
the thick solid line shows qEAM results and circles de- 
note QM results of Section II. We used a simulation cell 
containing 54 atoms for the perfect crystal case with pe- 
riodic boundary conditions. Vacancy formation enthalpy 
is defined in the same way as CyadP) [eq. (p^)] replacing 
energy with enthalpy. In Figure 6(b) we plot the va- 
cancy formation enthalpy hyac with respect to pressure. 
The difference between the vacancy formation enthalpy 
as obtained by qEAM FF and QM, is due to the differ- 
ence in the vacancy formation volume obtained from the 
two methods. The vacancy formation volume \vvac{P)] 
is, again, defined in the same way as evac{P) [eq. (M)] 
replacing energy with volume. In Figure 7 we plot the 
vacancy formation volume {fl^ vac), with respect to pres- 
sure. While the pressure dependence of Vvac{P) calcu- 
lated using the qEAM FF is very similar to the MB-PS 
calculation, our Force Field overestimates the vacancy 
formation volume resulting in higher vacancy formation 
enthalpy for compressed states. 

In order to calculate the vacancy formation enthalpy 
at finite temperatures as a function of pressure we per- 
formed NPT MD simulations using a cell containing 
N — 1458 atoms with periodic boundary conditions. We 
performed simulations at 9 different volumes (19, 18.36, 
18, 17, 16, 15, 14, 13, 11 A^); for each volume we started 
with the relaxed structure a T = K and heat the sys- 
tem in 100 K steps; for each temperature we performed 25 
ps MD simulations and used the last 20 ps for measure- 
ments. In Figure 8 we show vacancy formation enthalpy 
as a function of pressure for T — OK (full atomic relax- 
ation), T = 300 K, the volume- relaxed enthalpy is also 
shown for comparison; the zero pressure values are also 
shown in Table The fundamental data used to calcu- 
late vacancy energy, enthalpy pad volume can be found 
in the supplementary material.E£l 

A very important quantity, which determines vacancy 
mobility, is the vacancy migration energy barrier (i?™/)- 
We calculate E^J'^ using the qEAM FF by marching a 
nearest neighbor atom towards the position of the va- 
cancy in short steps. At each step the position of the 
marching atom is fixed, as well as that of a distant, ref- 
erence, atom, and the positions of all the other atoms are 
relaxed at constant pressure. In this way we obtain the 
optimum migration path and energy as a function of dis- 
placement. In Figure 9 we show the energy as a function 
the position of the marching atom for zero pressure and 
zero temperature; we obtain a vacancy migration energy 
^vac ~ 1-093 eV. The activation energy for self diffusion 
is defined as Q = Ei^^ + El^^s , using the qEAM FF we 
obtain Q = 4.028 eV in good agreement with the expe&i 
imental value 3.8 ± 0.3 eV,EJ and ah initio calculationsci 
which give 3.82 eV, see Table 0. At T=300 K we obtain 
a vacancy migration energy of 1.1 ± 0.5 eV, see Table 
IV, very similar to the zero temperature value. 
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C. Energetics of homogeneously sheared bcc crystal 

Zero temperature, first principles energetics of a homo- 
geneously sheared bcc crystal in the observed twinning 
mode was also included in the qEAM training set. 

The ideal shear strength is defined to be the stress 
separating elastic and plastic deformation when a homo- 
geneous shear is applied to a perfect crystal. It gives an 
upper bound for the shear strength of the material. The 
shear transformation is in the direction of the nbaei'ved 
twining mode and deforms the crystal into itselfES For 
bcc crystaLwe use the following transformation of the cell 
vectorsw'cj 



[111], 



b4[lll] + -^[Tll], 



c=^[lll], 



(15) 



when the shear variable s is equal to the twinning shear 
s = Stw = 2^^/^ the lattice vectors [a — 1/3[212], b — 
1/3[122] and c = ^[111]] form a bcc structure, twin of 
the initial one. 

In this way one can calculate the energy along the shear 
path. 



Wis) = e{V, s) - e{V, s = 0), 



(16) 



where e(V, s) is the energy per atom of the deformed 
system and e{V,s — 0) is the perfect crystal energy. 
The energy barrier associated with this transformation 
is Wmax = W{s = 0.5). The corresponding stress is de- 
fined as: 



r{s) 



1 dW{s) 
V ds 



(17) 



The ideal shear strength {Tmax) is defined as the maxi- 
mum stress along the path. In Figure 10 we show energy 
and stress as a function of shear using the qEAM FF for 
zero pressure volume, seCjalso Ref. ttq. 

Soderlind and Moriartya calculated W{s) and r(s) for 
Ta at different volumes, from first principles. In develop- 
ing the qEAM FF we used Wmax for ¥=17.61861^ and 
¥=10.9091^ as part of the training set. In Table |^ we 
show a comparison between the first principles resultsEJ 
and the ones obtained using the qEAM FF. We can see 
that the qEAM results are is in very good agreement with 
the ah initio calculations. 



D. Surface energy 

The unrelaxed (100) surface energy using the qEAM 
FF is 1.971 J/m^, lower than the first principles values 
of 3.097 J/m2 of Ref. M and 3.23 J/m of Ref. IM Low 



fields, see Ref. |l|. The zero temperature experimental 
estimate of the surfacej-energy (averaged over different 
surfaces) is 2.902 J/m^.E3 



IV. THERMAL EXPANSION 

Thermal expansivity is an important materials prop- 
erty that can be calculated directly from MD simulations. 
In order to calculate the lattice constant as a function of 
temperature for zero pressure we performed NPT MD 
simulations with a computational cell containing 1024 
atoms increasing the temperature by 100 K every 25 ps 
at zero applied pressure. For each temperature the first 5 
ps were taken as thermalization and the remaining 20 ps 
were used for measurements.tS ah initio MD simulations 
are very time consuming and only feasible for small sys- 
tems and short times; thus in order to calculate the ther- 
mal expansion from first principles we take the following 
approach.L3 The Helmholtz free energy can be written 



F{V,T) ^ E^{V) + Fei{V,T) + F^,b, 



(18) 



where E^iV) is the zero temperature energy (Section II), 
Ff,i{y,T) is the electronic contribution and Fmi, is the 
vibrational part of the free energy. The electronic contri- 
bution is obtained using quantum statistical mechanics. 
The vibraticpal part is obtained using the particle in a 
cell method,Ej further detail of our calculations can be 
found in Ref. ^ This calculations were performed us- 
ing the mixed basis pseudopotential method and a cell 
containing 54 atoms. 

In Figure 11(a) we show the thermal expansivity 
(a(T) = \/VdV/dT) as a function of temperature ob- 
tained form our MD simulations (circles), first principles 
results (dashed-dotted line), a|S_.well as the experimental 
values (full and dashed lines). l3 Figure 11(b) shows the 
linear thermal expansion (a — ao)/ao (where oq is the 
T= 300 K lattice constant) as a function of tempera- 
ture obtained using tlie qEAM FF (circles) as well as the 
experimental values.c3 The high temperature experimen- 
tal results (shown dashed lines) represent provisional 
experimental data.c2l Both MD and the particle in cell 
methods are based on classical mechanics so none of the 
them are expected to capture the low temperature behav- 
ior where the differences between quantum and classical 
statistical mechanics are important. The force field re- 
sults agree with experimental data well; for example the 
qEAM FF underestimates the change in lattice parame- 
ter from T = 300 K to T = 2000 K by less than 0.2%. 
This is an important result since the thermal expansion 
in related to anharmonicities of the internal energy. Our 
ah initio data also agree well with the experimental re- 
sults; we slightly overestimate thermal expansivity for 
temperatures lower than T=2000 K. 



surface energy is a common problem in EAM-like force 
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V. MELTING CURVE OF Ta 

We studied melting of Ta using one phase and two- 
phase MD simulations with the qEAM force field. 

In Figure 12 we show the results of zero pressure heat- 
ing a solid sample until it melts and then cooling the melt. 
The system contains 1024 atoms with periodic boundary 
conditions with a heating and cooling rate of 100 K per 
25 ps. In Figure 12 we plot enthalpy [Figure 12(a)] and 
volume [Figure 12(b)] as a function of temperature. The 
solid superheats on heating and the liquid undercook as 
expected for a small, periodic system with no defects to 
act as nucleation centers and high heating and cooling 
rates. 

In order to overcome some of these problems we cal- 
culate the melting temperature using the "two phase 
technique". We place pre-equilibrated liquid and pre- 
equilibrated solid samples in a single computational cell. 
Once this initial configuration is built we perform TPN 
MD (using a Hoover thermostat and Rahman-Parinello 
barostat) simulation and observe which phase grows. If 
the simulation temperature T is lower than the melting 
temperature at the given pressure [TmeitiP)] the liquid 
will start crystallizing, on the other hand the solid will 
meh if T> r™e/t(/'). 

Given that the liquid-solid phase transition is first or- 
der, i.e. the enthalpy of the solid and liquid in equilib- 
rium differ by a finite amount, it is very easy to know 
whether the system is crystallizing or melting by analyz- 
ing the time evolution of the total potential energy during 
the MD run. In Figure 13 we show potential energy as a 
function of time for zero pressure two-phase simulations 
at different temperatures. We can see that at T=3100 K 
(dotted line in Figure 13) the potential energy decreases 
with time; this means that the solid phase is growing 
and Tmeit > 3100 K. On the other hand, for T = 3200 
K the energy grows, the system is, then, melting and 
Tmeit < 3200 K. At T = 3150 K the energy is rather 
constant and both phases are close to equilibrium. This 
value for the zero pressure melting is in very good agree- 
ment with pcperimental results which range from 3213 K 
to 3273 K.cd This is a very important validation of the 
FF, taking into account that only zero temperature data 
was used in its development. 

The slope of the melting curve is given by the Clausius 
Clapeyron equation: 



dP 
'dT 



1 AH 

T Ay ' 



(19) 



where AH and AV are the enthalpy and volume differ- 
ence between the liquid and solid in equilibrium respec- 
tively. From our MD runs, see Figure 12, the slope of the 
melting curve at zero pressure is dT/dP = 92.8 K/GPa, 
larger than the the experimental value dT/dP = 60 ± 10 
K/GPa.y 

Using the two phase simulation procedure we calcu- 
lated the melting temperature for various pressures up 



to 300 GPa. In Figure 14 we show the melting curve for 
Ta, see also Ref. To the best of our knowledge this 
is the first calculation of melting temperature in Ta for 
such a— wide pressure range. Zero pressure experimental 
valuealJ are also shown in Figure 14 as empty diamonds. 
High pressure melting of Ta ha|S-|been studied experi- 
mentally via shock compression;E2l melting is identified 
as a change in the velocity of the rarefication wave (from 
the longitudinal to the bulk sound velocity). The transi- 
tion was found to^occur in the pressure range from ^^250 
GPa to 295 GPa.E3 The calculation of the temperature 
in shock experiments (i.e. along the Hugoniot) is diffi- 
cult; the electronic contribution to the specific beat has a 
very strong effect on the melting temperature.c3 Simple 
models for the electronic behavior lead to very differ- 
ent temperatures: using the free electron gas model the 
melting temperature is ^10000 K while considering band 
electrons give Tmeit ^ 7500 K, see Ref. These exper- 
imental values are shown in Figure 14 are empty circles. 
Using the accurate ab initio thermal equation of state 
obtained using the methods described above and the 
Rankine-Hugoniot equation Cohen and Giilseren calcu- 
lated pressure-volume and teniperature-pressure curves 
for Ta under shock conditions.L3 This calculation leads 
to a melting temperature (using the experimental melt- 
ing pressure Pme;t=375 GPa) of 8150 K (see square in 
Figure 14). Our MD results are in good agreement with 
the high pressure calculation of the melting temperature 
based on shock experiments and ab initio calculations. 



VI. CONCLUSIONS 

We have presented here a general strategy to derive 
accurate classical force fields based on ab initio QM me- 
chanical calculations for metallic systems. Force fields al- 
low calculations on systems containing millions of atoms, 
providing a means to study from an atomistic point of 
view a variety of processes relevant to the mechanical 
and thermodynamical properties of metals, such as phase 
transitions, dislocations dynamics and interactions, fail- 
ure, crack propagation, etc. 

We showed that the qEAM FF describes with good 
accuracy the EOS of Ta for different crystal phases [bcc 
(coordination number 8), fee (coordination number 12), 
and A15 (mixed coordination numbers)], and also elas- 
tic constants, vacancy formation energy, and energetics 
of the deformed bcc lattice in the twinning direction. A 
critical point in our approach is that in developing the 
force field we use not only the thermodynamically sta- 
ble phase but also high energy phases and large strains; 
a force field that can correctly describe such structures 
should be appropriate for simulations of defects and non- 
equilibrium processes. 

The large amount of QM data used to derive our force 
field gives an important measure of its quality. 

We used the qEAM FF with MD to calculate the melt- 
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ing curve for Ta in a wide pressure range. The zero pres- 
sure melting temperature obtained from our simulations 
Tmeit = 3150 ± 50 K is in very good agreement with the 
experimental result of 3290± 100 K. This is an important 
validation of our approach given the fact that the qEAM 
FF was derived using zero temperature data. 

First-priciples-based force fields represent an impor- 
tant step in ab initio multiscale modeling of materials. 
We have used the qEAM FF to study spall failure ,ca 
crack propagation, and dislocation properties such as 
core structure gmd energy, Peierls stress, and kink for- 
mation energiesE^ with systems containing as many as 
50,000 atoms. We used such calculations with a mi- 
cromechanical mcw^l of plasticity developed by Stainer, 
Cuitifio and Ortialj to develop a multiscale model of sin- 
gle crystal plasticity in Ta.E3 We are currently performing 
shock simulations in systems containing more than a mil- 
lion atoms using a parallel MD code; for a system con- 
taining 1,098,500 atoms 1 ps of MD simulations (1000 
steps) takes approximately one hour on an SGI Origin 
2000 computer using 128 250-MHz RIOOOO processors. 
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TABLE I. EOS parameters for bcc Tantalum. 





Vo{A^) 


BT(GPa) 




cii (GPa) 


Cl2 (GPi 

159.8 
137.57 
163 

:) 

IK) 
160.94 




Theory (0 K) 


This work (LAPW-GGA) 
This work (qEAM EE) 
EP LMTO GGA SC ^ 


18.33 
18.36 
17.68 


188.27 
183.04 
203 


4.08 
4.16 


245.18 
272.54 
281 




Theory (300 ¥ 


This work (qEAM EE) 


18.4 


176 


4.9 






Experiment (30( 


Diamond Anvil Cell " 
Ultrasonic 


18.04 


194.7±4.8 


3.4 


266 



^EuU potential linear mufli*i-tin orbital calculations, P. 
Soderlind and J. A. M|«riarty.tl 
'^H. Cynn and C. YooEJ 

'^Adiabaticieiastic constants at 25° C, Katahara, Manghnani, 
and Fisher O 



TABLE II. qEAM paramters for Ta: two body term. The 
units of bi are eVA"<^+^'. 





6o 




62 


fe3 


4.81253968 


6.50281587 


-11.26455130 


8.01451544 


-2.97299223 




bs 




fey 




0.60004206 


-0.06222106 


0.00258801 


-0.00000504 





TABLE III. qEAM paramters for Ta: embedding energy 



ai 



02 



(1/A) 



/3 



(A) 



0.07293238 



0.15781672 



21.79609053 



7.79329426 



3.32389219 



E,oh(eSI) 



Bt (GPa 



A 



h 



8.15420437 



183.035404 



0.20782789 



-0.00717801 



-0.00000504 



TABLE IV. Volume-relaxed (e^ac), full-relaxed Vacancy 
formation {e'^^), vacancy migration energies (e^J*), and ac- 
tivation energy for self diffusion (Q). 





^vac (^^) 


eSL(eV) 


Cac^(eV) 


Q (eV) 




Theory (0 K) | 


This work (qEAM EE) 


3.10 


2.935 


1.093 


4.028 


This work (MB-PS) 


3.25 








EP LMTO GGA SC 


3.74 


2.20 






Plane waves LDA 


3.51 


2.99 


0.83 


3.82 


EP LMTO LDA 


3.49 










Theory and Experiment (300 K) | 


This work (qEAM EE) 
Experiment 




3.0 ± 0.05 


1.1 ± 0.05 


4.1 ± 0.1 




2.8 ± 0.6 




3.8 ± 0.3 



''P. Soderlind and J. A. Moriartjj 
''Satta, Willaime, and Gironcoli.l 
'^Unrelaxed value by Korhonen, Puska, and Nieminer.Ej 
"^Ref. ^ for vacancy energy and Ref. ^ for activation energy 
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TABLE V. Shear deformation in the observed twinning 
mode in Ta. 



Volume (A^) 


Wr^a. (eV) 


'Tmax 


Wr^a. (eV) 






This work (qEAM FF) 


FP LMTO GGA SC ^ 


18.36 


0.188 


7.14 






17.618602 


0.2 


8.0 


0.194 


7.37 


15.143996 


0.26 


12.05 


0.276 


12.4 


10.9090116 


0.43 


28.2 


0.566 


36.2 



""P. Soderlind and J. A. MoriartyJ 



FIGURE CAPTIONS 

Figure 1. Zero temperature EOS for bee Ta, LAPW 
GGA and qEAM FF results. Energy [Figure 1(a)] and 
pressure [Figure 1(b)] as a function of volume. Circles 
denote LAPW GGA results and lines show qEAM FF 
results. 

Figure 2. Zero temperature EOS for fee Ta, LAPW 
GGA and qEAM FF results. Energy [Figure 2(a)] and 
pressure [Figure 2(b)] as a function of volume. Circles 
denote LAPW GGA results and lines show qEAM FF 
results. The inset of Figure 2(a) shows fcc-bcc energy 
difference as a function of volume (circles joined by dots 
denote LAPW GGA results and line qEAM FF). 

Figure 3. Zero temperature EOS for A15 Ta, ab initio 
and qEAM FF results. Energy [Figure 3(a)] and pres- 
sure [Figure 3(b)] as a function of volume. Open circles 
with dotted line denote full potential muffin-tin orbital 
calculations results Ref. |^; filled. circles show the sum of 
Al5-bcc energy difference fromld and our calculation of 
bcc energy using LAPW GGA method (section II). The 
full line shows qEAM FF results. The inset of Figure 3(a) 
shows A15-bcc energy difference as a function of volume 
(circles with dotted line denote QM resultil and the line 
qEAM FF). 

Figure 4. Zero temperature EOS for hep Ta, LAPW 
GGA and qEAM FF results. Energy [Figure 4(a)] and 
pressure [Figure 4(b)] as a function of volume. Circles 
denote LAPW GGA results, section II, and lines show 
qEAM FF results. The inset of Figure 4(a) shows hcp- 
bcc energy difference as a function of volume (circles with 
dotted line denote LAPW GGA results and line qEAM 
FF). 

Figure 5. Zero temperature elastic constants for Ta, 
LAPW GGA and qEAM FF results. Circles show bulk 
modulus [(cii -f 2ci2)/3]; diamonds show C44 and squares 
represent Cs — (cn — ci2)/2. qEAM FF results are shown 
with filled symbols and full lines and ab initio LAPW 
results with open symbols and dashed lines. 

Figure 6. Volume relaxed vacancy formation energy 

(a) and enthalpy (b) as function of pressure. Dashed 
lines represent ab initio MB-PS results (section II) and 
full lines show qEAM results. 

Figure 7. Vacancy formation volume as function of 
pressure. Dashed lines represent ab initio MB-PS results 
(section II) and full lines show qEAM results. 

Figure 8. Vacancy formation enthalpy as function of 
pressure using the qEAM FF. Solid line shows the volume 
relaxed result; the dotted line the T=0 K fully relaxed 
results; dashed line is T=300 K result. 

Figure 9. Vacancy migration energy using qEAM FF. 
Energy as a function of position of the marching atom at 
T=0 K and zero pressure. The vacancy migration energy 
is 1.093 eV. 

Figure 10. Ideal shear strength of Ta using qEAM FF 
at zero temperature and volume V=18.36 A^. We show 
energy W{s) [Figure 10 (a)] and stress t(s) [Figure 10 

(b) ] as a function of shear. 
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Figure 11. Thermal expansion in Ta. (a) thermal 
expansivity as a function of temperature; circles repre- 
sent qEAM FF results, the dashed-dotted lines shows 
mixed basis pseudopotantial calculations using the par- 
ticle in a cell method, and the line denotes experimental 
results from Ref. (b) Linear thermal expansion of 
Ta. (a — ao)/ao as a function of temperature; circles rep- 
resents qEAM FF results and line denote experimental 
results from Ref. |4^. The high temperature^xperimental 
data (dashed lines) are provisional values.cHi 

Figure 12. Melting of tantalum using the qEAM FF. 
Enthalpy (a) and volume (b) as a function of tempera- 
ture for zero pressure; heating of bcc Ta (lower branches) 
and cooling of liquid Ta (higher branches). Heating and 
cooling rates are 100 K per 25 ps. 

Figure 13. Ta melting using the qEAM FF. Two phase 
simulations. Time evolution of the potential energy in 
TPN MD at zero pressure for different temperatures. For 
T=3100 K (dotted line) the potential energy decreases 
with time, i.e. the system is crystallizing; for T=3200 K 
(dashed line) the potential energy grows denoting melt- 
ing. For T = 3150 ~ Tmeit{P = 0) the potential energy 
remains constant; in this case the solid and liquid phases 
coexist in equilibrium. 

Figure 14. Ta Melting qEAM FF. Melting curve for 
Ta up to P= 300 GPa obtained using the two-phase sim- 
ulation technique. We also show the experimental zero 
pressure melting ptemperatureO (circles) and the results 
of shock melting.E3 
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